function EL = find_lines (X,Y,Z)

A = [X(:) Y(:) Z(:)];
A(end+1,:) = [X(1) Y(1) Z(1)];
rest_inds = 1:length(X);

DX = diff(A(:,1));
DY = diff(A(:,2));
DZ = diff(A(:,3));

D = (DX.^2+DY.^2+DZ.^2).^0.5; % Distance

el_count = 1;

while(length(rest_inds) > 1)
    
    rest_points = A(rest_inds,:);
    spoint = rest_points(1,:);
    dtp =  distancePoint(spoint,rest_points(2:end,:));
    
    if(min(dtp) > 10)
        rest_inds(1) = [];
        continue;
    end
    
    [dtps ind] = sort(dtp);
    
    start_inds = [rest_inds(1) rest_inds(ind(1)+1)];
    rest_inds([1 ind(1)+1]) = [];
%     rest_inds = rest_inds([1:ind(1)-1 (ind(1)+2):length(rest_inds)]);
    b = length(rest_inds);
    for i=1:b
        
        points = A(start_inds,:);
        rest_points = A(rest_inds,:);
        
        dtl = distanceToLine(points,rest_points);
        
        if(min(dtl) > 3)
            break;
        end
        
        [ds indd] = sort(dtl);
        
        for j=1:length(indd)
            
            dtp =  distancePoint(rest_points(indd(j),:),points);

            
            [dtps] = sort(dtp);
            
            if(dtps(1) < 10 && ds(j) < 3 && dtps(1) > 4)
                
                start_inds(end+1) = rest_inds(indd(j));
                rest_inds(indd(j)) = [];
                
                break;
            end
            
        end
        
        
    end
    
    if(length(start_inds) >=5 )
        
        EL(el_count) = {start_inds};
        el_count = el_count + 1;
        
%     elseif(length(rest_inds)>5)
%         rest_inds = [rest_inds start_inds];
    end
    
    
end

end




function d = distanceToLine(points,rest_points)

[coeff] = princomp(points);
meanX = mean(points,1);
dirVect = coeff(:,1);
meanX2points = (rest_points-repmat(meanX,size(rest_points,1),1));
D = dot(meanX2points',repmat(dirVect',size(rest_points,1),1)');

dVector = (dirVect*D)'-meanX2points;
% D = D + rest_points-repmat(meanX,size(rest_points,1),1);

d = (dVector(:,1).^2+dVector(:,2).^2+dVector(:,3).^2).^0.5;

end

function d = distancePoint(p1,x0)

D = x0 - repmat(p1,size(x0,1),1);

d = (D(:,1).^2+D(:,2).^2+D(:,3).^2).^0.5;

end


